function [ Qn, Grad_Qn, Var_GMM] = GMM_Fn_Mar2022(starts,Cons,G_ij,y_in)

    G_ji = Cons.Flip_ij*G_ij;
    G_ij_doti = Cons.Trans_doti*G_ij;
    G_ji_doti = Cons.Trans_doti*G_ji;

    obs = size(Cons.X_i,1);
    
    dim_BGD = 4*Cons.dimX+3;
    BGD_starts = starts(1:dim_BGD);
    A_starts = starts(4+4*Cons.dimX:dim_BGD+Cons.n);

    %ba = BGD_starts(1);
    gamma = BGD_starts(2:2+Cons.dimX);
    delta = BGD_starts(3+Cons.dimX:dim_BGD);
    A_i = Cons.I_i*A_starts;
    A_j = Cons.I_j*A_starts;
    A_j_doti = Cons.Trans_doti*A_j;
   
    
    RHS_doti = [G_ji_doti, Cons.X_j_doti, A_j_doti, bsxfun(@times,Cons.X_i,Cons.X_j_doti), bsxfun(@times,Cons.X_i,A_j_doti), bsxfun(@times,A_i,Cons.X_j_doti), bsxfun(@times,A_i,A_j_doti)];
    C_hat_doti = -G_ij_doti + RHS_doti*BGD_starts;
    
    %% First Moment Conditions (fixed effects)
    %Z1 = Cons.I_j;       
    MomCon1 = bsxfun(@times,C_hat_doti,Cons.I_j);
    
    % Gradient
    %D_A_j = Cons.I_j;
    D_A_j_doti = Cons.Trans_doti*Cons.I_j;
    
    dC_hat_doti = [RHS_doti, bsxfun(@times,(gamma(1+Cons.dimX)*ones(size(Cons.I_j,1),1) + Cons.X_i*delta(1+Cons.dimX:2*Cons.dimX)+A_i*delta(3*Cons.dimX+1)),D_A_j_doti) + bsxfun(@times,(Cons.X_j_doti*delta(2*Cons.dimX+1:3*Cons.dimX) + A_j_doti*delta(3*Cons.dimX+1)),Cons.I_i)];
    
    d_MomCon1 = sparse(Cons.I_j'*dC_hat_doti);
    
    %% Second Moment Conditions (independence and variance of A)
    Z2 = [ones(obs,1), Cons.X_i, Cons.X_j, bsxfun(@times,Cons.X_i,Cons.X_j)];
    MomCon2 = [bsxfun(@times,A_j,Z2), bsxfun(@times,bsxfun(@times,A_j,A_i),Z2), bsxfun(@times,A_j.^2 - 1,Z2)];
        
    % Gradient
    d_MomCon2 = [Z2'*Cons.I_j; (bsxfun(@times,Z2,A_i)'*Cons.I_j + bsxfun(@times,Z2,A_j)'*Cons.I_i); 2*bsxfun(@times,Z2,A_j)'*Cons.I_j];
    d_MomCon2 = [sparse(size(d_MomCon2,1),dim_BGD), d_MomCon2];     
    
    
    %% Third Moment Conditions (Independence of X and A from C)       
    % Instruments for third set of moments
    sum_A = sum(bsxfun(@times,A_starts',Cons.BlockOnes),2);
    mean_A_k_less_ij = bsxfun(@rdivide,sum_A - A_i - A_j,Cons.I_i*Cons.I_sum1');

    Z3 = sparse([RHS_doti(:,2:size(RHS_doti,2)), Cons.X_j_mean_Xk_less_ij, Cons.X_i_mean_Xk_less_ij, bsxfun(@times,A_j,mean_A_k_less_ij), bsxfun(@times,A_i,mean_A_k_less_ij), bsxfun(@times,Cons.X_j,mean_A_k_less_ij), bsxfun(@times,Cons.X_i,mean_A_k_less_ij), bsxfun(@times,A_j,Cons.mean_Xk_less_ij), bsxfun(@times,A_i,Cons.mean_Xk_less_ij)]);

    MomCon3 = bsxfun(@times, C_hat_doti, Z3);

    % Gradient
     d_mean_A_k_less_ij = bsxfun(@rdivide,Cons.BlockOnes - Cons.I_i - Cons.I_j,Cons.I_sum1);

     dZ3_Chat_doti = [sparse(Cons.dimX,Cons.n); ...
         C_hat_doti'*D_A_j_doti;
         sparse(Cons.dimX,Cons.n); ...
         bsxfun(@times,C_hat_doti,Cons.X_i)'*D_A_j_doti; ...
         bsxfun(@times,C_hat_doti,Cons.X_j_doti)'*Cons.I_i; ...
         bsxfun(@times,C_hat_doti,A_j_doti)'*Cons.I_i + bsxfun(@times,C_hat_doti,A_i)'*D_A_j_doti; ...
         sparse(2*Cons.dimX,Cons.n); ...
         C_hat_doti'*(bsxfun(@times,A_j,d_mean_A_k_less_ij) + bsxfun(@times,mean_A_k_less_ij,Cons.I_j) ); ...
         C_hat_doti'*(bsxfun(@times,A_i,d_mean_A_k_less_ij) + bsxfun(@times,mean_A_k_less_ij,Cons.I_i) ); ...
         bsxfun(@times,C_hat_doti,Cons.X_j)'*d_mean_A_k_less_ij; ...
         bsxfun(@times,C_hat_doti,Cons.X_i)'*d_mean_A_k_less_ij; ...
         bsxfun(@times,C_hat_doti,Cons.mean_Xk_less_ij)'*Cons.I_j; ...
         bsxfun(@times,C_hat_doti,Cons.mean_Xk_less_ij)'*Cons.I_i];
     
     
     dC_hat_doti_Z3 = Z3'*dC_hat_doti;
   
     d_MomCon3 = [sparse(size(dZ3_Chat_doti,1),dim_BGD), dZ3_Chat_doti] + dC_hat_doti_Z3;    

    %% Conditions and Derivatives
    G1 = [MomCon3, MomCon1, MomCon2];
    Grad1 = [d_MomCon3; d_MomCon1; d_MomCon2]/obs;
    
    W = speye(size(G1,2));
    G1_mean = mean(G1);
    Qn = full(G1_mean*W*G1_mean');
    Grad_Qn = full(2*Grad1'*W*G1_mean');    
    
    %%% Calculate Variance if called for
    if nargout>2
        W_ij = exp(G_ij) + exp(G_ji);
        W_ij = Cons.W_trans*W_ij;
        W_ij_rep = zeros(Cons.n,Cons.n);
        k = 0;
        for s=1:Cons.schools
            W_ij_rep(k+1:k+Cons.sch_size(s),k+1:k+Cons.sch_size(s)) = reshape(W_ij(k+1:k+Cons.sch_size(s)^2,1),Cons.sch_size(s),Cons.sch_size(s));
            k = k + Cons.sch_size(s);
        end    
        W_ij_MAT = bsxfun(@rdivide,W_ij_rep,sum(W_ij_rep,2));
        
        Trans_mat = (obs/Cons.n)*Cons.I_i./sqrt((repmat(Cons.I_sum,obs,1)));
        
        d_A = Cons.I_i;
        d_W_ij_A = Cons.I_i*W_ij_MAT;
        d_W2_ij_A = Cons.I_i*W_ij_MAT^2;
        
        dim_alpha = 60;
        if size(starts,1)==dim_BGD+Cons.n+2*dim_alpha
            dim_beta = 0;
        else
            dim_beta = 24;
        end
        
        %% Education
        alpha_educ = starts(dim_BGD+Cons.n+1:dim_BGD+Cons.n+dim_alpha);
        alpha_educ_m1 = alpha_educ(1:9);
        dim_X_m1 = size(alpha_educ_m1,1);
        alpha_educ_m2 = alpha_educ(10:24);
        dim_X_m2 = size(alpha_educ_m2,1);
        alpha_educ_m3 = alpha_educ(25:39);
        dim_X_m3 = size(alpha_educ_m3,1);
        alpha_educ_m4 = alpha_educ(40:60);
        dim_X_m4 = size(alpha_educ_m4,1);         
  
        LMH_educ = [Cons.L_educ, Cons.M_educ, Cons.H_educ];
        
        X_0 = [ones(Cons.n,1), Cons.P, W_ij_MAT*Cons.P];
        X_0_educ_LMH = [];
        for z=1:3
            X_0_educ_LMH = [X_0_educ_LMH, bsxfun(@times,LMH_educ, X_0(:,z))];
        end
        
        % M1
        X_educ_m1 = Trans_mat*X_0_educ_LMH;
        Ehat_educ_m1= Trans_mat*y_in.y1_educ_m1 - X_educ_m1*alpha_educ_m1;
        MC_educ_m1 = bsxfun(@times,Ehat_educ_m1,X_educ_m1);

        dEhat_educ_m1_alpha = -X_educ_m1;
        dEhat_educ_m1 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), dEhat_educ_m1_alpha, sparse(obs,dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_beta)];
        X_dEhat_educ_m1 = X_educ_m1'*dEhat_educ_m1;
        
        d_educ_m1 = X_dEhat_educ_m1;
        
        % M2
        X_educ_m2 = Trans_mat*[X_0_educ_LMH, bsxfun(@times,LMH_educ,A_starts), bsxfun(@times,LMH_educ,W_ij_MAT*A_starts)];
        Ehat_educ_m2= Trans_mat*y_in.y1_educ_m2 - X_educ_m2*alpha_educ_m2;
        MC_educ_m2 = bsxfun(@times,Ehat_educ_m2,X_educ_m2);
        
        dEhat_educ_m2_alpha = -X_educ_m2;
        dEhat_educ_m2_A = -(bsxfun(@times,Trans_mat*LMH_educ(:,1),alpha_educ_m2(10)*d_A + alpha_educ_m2(13)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,2),alpha_educ_m2(11)*d_A + alpha_educ_m2(14)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,3),alpha_educ_m2(12)*d_A + alpha_educ_m2(15)*d_W_ij_A) );
        dEhat_educ_m2 = [sparse(obs,dim_BGD), dEhat_educ_m2_A, sparse(obs,dim_X_m1), dEhat_educ_m2_alpha, sparse(obs,dim_X_m3+dim_X_m4), sparse(obs,dim_beta)];
        X_dEhat_educ_m2 = X_educ_m2'*dEhat_educ_m2;
       
        dX_Ehat_educ_m2_A = [sparse(9,Cons.n); ...
            Ehat_educ_m2'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_A); Ehat_educ_m2'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_A); Ehat_educ_m2'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_A); ...
            Ehat_educ_m2'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_W_ij_A); Ehat_educ_m2'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_W_ij_A); ...
            Ehat_educ_m2'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_W_ij_A)];
        
        dX_Ehat_educ_m2 = [sparse(dim_X_m2,dim_BGD), dX_Ehat_educ_m2_A, sparse(dim_X_m2,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        
        d_educ_m2 = X_dEhat_educ_m2 + dX_Ehat_educ_m2;
                
        % M3
        X_educ_m3 = Trans_mat*[X_0_educ_LMH, bsxfun(@times,LMH_educ,Cons.y0_educ), bsxfun(@times,LMH_educ,W_ij_MAT*Cons.y0_educ)];
        Ehat_educ_m3 = Trans_mat*y_in.y1_educ_m3 - X_educ_m3*alpha_educ_m3;
        MC_educ_m3 = bsxfun(@times,Ehat_educ_m3,X_educ_m3);       
        
        dEhat_educ_m3_alpha = -X_educ_m3;
        dEhat_educ_m3 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), sparse(obs,dim_X_m1+dim_X_m2), dEhat_educ_m3_alpha, sparse(obs,dim_X_m4), sparse(obs,dim_beta)];
        X_dEhat_educ_m3 = X_educ_m3'*dEhat_educ_m3;
        
        d_educ_m3 = X_dEhat_educ_m3;
        
        % M4
        X_educ_m4 = Trans_mat*[X_0_educ_LMH, bsxfun(@times,LMH_educ,A_starts), bsxfun(@times,LMH_educ,W_ij_MAT*A_starts), bsxfun(@times,LMH_educ,Cons.y0_educ), bsxfun(@times,LMH_educ,W_ij_MAT*Cons.y0_educ)];
        Ehat_educ_m4= Trans_mat*y_in.y1_educ_m4 - X_educ_m4*alpha_educ_m4;
        MC_educ_m4 = bsxfun(@times,Ehat_educ_m4,X_educ_m4);
         
        dEhat_educ_m4_alpha = -X_educ_m4;
        dEhat_educ_m4_A = -(bsxfun(@times,Trans_mat*LMH_educ(:,1),alpha_educ_m4(10)*d_A + alpha_educ_m4(13)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,2),alpha_educ_m4(11)*d_A + alpha_educ_m4(14)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,3),alpha_educ_m4(12)*d_A + alpha_educ_m4(15)*d_W_ij_A) );
        dEhat_educ_m4 = [sparse(obs,dim_BGD), dEhat_educ_m4_A, sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3), dEhat_educ_m4_alpha, sparse(obs,dim_beta)];
        X_dEhat_educ_m4 = X_educ_m4'*dEhat_educ_m4;
       
        dX_Ehat_educ_m4_A = [sparse(9,Cons.n); Ehat_educ_m4'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_A); ...
            Ehat_educ_m4'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_A); Ehat_educ_m4'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_A); ...
            Ehat_educ_m4'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_W_ij_A); Ehat_educ_m4'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_W_ij_A); ...
            Ehat_educ_m4'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_W_ij_A); sparse(6,Cons.n)];
        dX_Ehat_educ_m4 = [sparse(dim_X_m4,dim_BGD), dX_Ehat_educ_m4_A, sparse(dim_X_m4,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        
        d_educ_m4 = X_dEhat_educ_m4 + dX_Ehat_educ_m4;
        
        
        
        %%% GENDER
        %% Gender
        alpha_gender = starts(dim_BGD+Cons.n+(dim_alpha+dim_beta)+1:dim_BGD+Cons.n+(dim_alpha+dim_beta)+dim_alpha);
        alpha_gender_m1 = alpha_gender(1:9);
        dim_X_m1 = size(alpha_gender_m1,1);
        alpha_gender_m2 = alpha_gender(10:24);
        dim_X_m2 = size(alpha_gender_m2,1);
        alpha_gender_m3 = alpha_gender(25:39);
        dim_X_m3 = size(alpha_gender_m3,1);
        alpha_gender_m4 = alpha_gender(40:60);
        dim_X_m4 = size(alpha_gender_m4,1);
        

        LMH_gender = [Cons.L_gender, Cons.M_gender, Cons.H_gender];

        X_0 = [ones(Cons.n,1), Cons.P, W_ij_MAT*Cons.P];
        X_0_gender_LMH = [];
        for z=1:3
            X_0_gender_LMH = [X_0_gender_LMH, bsxfun(@times,LMH_gender, X_0(:,z))];
        end
        
        % M1
        X_gender_m1 = Trans_mat*X_0_gender_LMH;
        Ehat_gender_m1= Trans_mat*y_in.y1_gender_m1 - X_gender_m1*alpha_gender_m1;
        MC_gender_m1 = bsxfun(@times,Ehat_gender_m1,X_gender_m1);
 
        dEhat_gender_m1_alpha = -X_gender_m1;
        dEhat_gender_m1 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), dEhat_gender_m1_alpha, sparse(obs,dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_beta)];
        X_dEhat_gender_m1 = X_gender_m1'*dEhat_gender_m1;
        
        d_gender_m1 = X_dEhat_gender_m1;
                

        % M2
        X_gender_m2 = Trans_mat*[X_0_gender_LMH, bsxfun(@times,LMH_gender,A_starts), bsxfun(@times,LMH_gender,W_ij_MAT*A_starts)];
        Ehat_gender_m2= Trans_mat*y_in.y1_gender_m2 - X_gender_m2*alpha_gender_m2;
        MC_gender_m2 = bsxfun(@times,Ehat_gender_m2,X_gender_m2);
        
        dEhat_gender_m2_alpha = -X_gender_m2;
        dEhat_gender_m2_A = -(bsxfun(@times,Trans_mat*LMH_gender(:,1),alpha_gender_m2(10)*d_A + alpha_gender_m2(13)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,2),alpha_gender_m2(11)*d_A + alpha_gender_m2(14)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,3),alpha_gender_m2(12)*d_A + alpha_gender_m2(15)*d_W_ij_A) );
        dEhat_gender_m2 = [sparse(obs,dim_BGD), dEhat_gender_m2_A, sparse(obs,dim_X_m1), dEhat_gender_m2_alpha, sparse(obs,dim_X_m3+dim_X_m4), sparse(obs,dim_beta)];
        X_dEhat_gender_m2 = X_gender_m2'*dEhat_gender_m2;
       
        dX_Ehat_gender_m2_A = [sparse(9,Cons.n); ...
            Ehat_gender_m2'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_A); Ehat_gender_m2'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_A); Ehat_gender_m2'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_A); ...
            Ehat_gender_m2'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_W_ij_A); Ehat_gender_m2'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_W_ij_A); ...
            Ehat_gender_m2'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_W_ij_A)];
        
        dX_Ehat_gender_m2 = [sparse(dim_X_m2,dim_BGD), dX_Ehat_gender_m2_A, sparse(dim_X_m2,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        
        d_gender_m2 = X_dEhat_gender_m2 + dX_Ehat_gender_m2;

        % M3
        X_gender_m3 = Trans_mat*[X_0_gender_LMH, bsxfun(@times,LMH_gender,Cons.y0_gender), bsxfun(@times,LMH_gender,W_ij_MAT*Cons.y0_gender)];
        Ehat_gender_m3 = Trans_mat*y_in.y1_gender_m3 - X_gender_m3*alpha_gender_m3;
        MC_gender_m3 = bsxfun(@times,Ehat_gender_m3,X_gender_m3);       
        
        dEhat_gender_m3_alpha = -X_gender_m3;
        dEhat_gender_m3 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), sparse(obs,dim_X_m1+dim_X_m2), dEhat_gender_m3_alpha, sparse(obs,dim_X_m4), sparse(obs,dim_beta)];
        X_dEhat_gender_m3 = X_gender_m3'*dEhat_gender_m3;
        
        d_gender_m3 = X_dEhat_gender_m3;
        
        % M4
        X_gender_m4 = Trans_mat*[X_0_gender_LMH, bsxfun(@times,LMH_gender,A_starts), bsxfun(@times,LMH_gender,W_ij_MAT*A_starts), bsxfun(@times,LMH_gender,Cons.y0_gender), bsxfun(@times,LMH_gender,W_ij_MAT*Cons.y0_gender)];
        Ehat_gender_m4= Trans_mat*y_in.y1_gender_m4 - X_gender_m4*alpha_gender_m4;
        MC_gender_m4 = bsxfun(@times,Ehat_gender_m4,X_gender_m4);
         
        dEhat_gender_m4_alpha = -X_gender_m4;
        dEhat_gender_m4_A = -(bsxfun(@times,Trans_mat*LMH_gender(:,1),alpha_gender_m4(10)*d_A + alpha_gender_m4(13)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,2),alpha_gender_m4(11)*d_A + alpha_gender_m4(14)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,3),alpha_gender_m4(12)*d_A + alpha_gender_m4(15)*d_W_ij_A) );
        dEhat_gender_m4 = [sparse(obs,dim_BGD), dEhat_gender_m4_A, sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3), dEhat_gender_m4_alpha, sparse(obs,dim_beta)];
        X_dEhat_gender_m4 = X_gender_m4'*dEhat_gender_m4;
       
        dX_Ehat_gender_m4_A = [sparse(9,Cons.n); Ehat_gender_m4'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_A); ...
            Ehat_gender_m4'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_A); Ehat_gender_m4'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_A); ...
            Ehat_gender_m4'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_W_ij_A); Ehat_gender_m4'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_W_ij_A); ...
            Ehat_gender_m4'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_W_ij_A); sparse(6,Cons.n)];
        
        dX_Ehat_gender_m4 = [sparse(dim_X_m4,dim_BGD), dX_Ehat_gender_m4_A, sparse(dim_X_m4,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        
        d_gender_m4 = X_dEhat_gender_m4 + dX_Ehat_gender_m4;        
        
        
        
        
        
        %%% 2SLS if called for
        if size(starts,1)>dim_BGD+Cons.n+2*dim_alpha
            
            %%% EDUC
            beta_educ = starts(dim_BGD+Cons.n+dim_alpha+1:dim_BGD+Cons.n+dim_alpha+dim_beta);
            dim_X_p1 = 0;
            dim_X_p2 = 0;
            %beta_educ_p3 = beta_educ(1:18);
            %dim_X_p3 = size(beta_educ_p3,1);
            dim_X_p3 = 0;
            beta_educ_p4 = beta_educ(1:24);
            dim_X_p4 = size(beta_educ_p4,1); 
                
            %         % P1
            %         X_educ_p1 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT*y_in.y1_educ_p1), X_educ_m1];
            %         Z_educ_p1 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*Cons.P), X_educ_m1];
            %         dim_Z_p1 = size(Z_educ_p1,2);
            %         Ehat_educ_p1= Trans_mat*y_in.y1_educ_p1 - X_educ_p1*beta_educ_p1;
            %         MC_educ_p1 = bsxfun(@times,Ehat_educ_p1,Z_educ_p1);
            %         
            %         dEhat_educ_p1_beta = -X_educ_p1;
            %         dEhat_educ_p1 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), dEhat_educ_p1_beta,  sparse(obs,dim_X_p2+dim_X_p3+dim_X_p4)];
            %         Z_dEhat_educ_p1 = Z_educ_p1'*dEhat_educ_p1;
            %         
            %         d_educ_p1 = Z_dEhat_educ_p1;
        

            %         % P2
            %         X_educ_p2 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT*y_in.y1_educ_p2), X_educ_m2];
            %         Z_educ_p2 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*Cons.P), Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*A_starts), X_educ_m2];
            %         dim_Z_p2 = size(Z_educ_p2,2);
            %         
            %         Ehat_educ_p2 = Trans_mat*y_in.y1_educ_p2 - X_educ_p2*beta_educ_p2;
            %         MC_educ_p2 = bsxfun(@times,Ehat_educ_p2,Z_educ_p2);
            %         
            %         dEhat_educ_p2_beta = -X_educ_p2;
            %         dEhat_educ_p2_A = -(bsxfun(@times,Trans_mat*LMH_educ(:,1),beta_educ_p2(13)*d_A + beta_educ_p2(16)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,2),beta_educ_p2(14)*d_A + beta_educ_p2(17)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,3),beta_educ_p2(15)*d_A + beta_educ_p2(18)*d_W_ij_A) );
            %         dEhat_educ_p2 = [sparse(obs,dim_BGD), dEhat_educ_p2_A, sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_X_p1), dEhat_educ_p2_beta, sparse(obs,dim_X_p3+dim_X_p4)];
            %         Z_dEhat_educ_p2 = Z_educ_p2'*dEhat_educ_p2;
            %         
            %         dZ_Ehat_educ_p2_A = [sparse(3,Cons.n); ...
            %             Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_W2_ij_A); Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_W2_ij_A); Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_W2_ij_A); ...
            %             sparse(9,Cons.n); ...
            %             Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_A); Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_A); Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_A); ...
            %             Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_W_ij_A); Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_W_ij_A); Ehat_educ_p2'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_W_ij_A)];
            % 
            %         dZ_Ehat_educ_p2 = [sparse(dim_Z_p2,dim_BGD), dZ_Ehat_educ_p2_A, sparse(dim_Z_p2,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
            %          
            %         d_educ_p2 = Z_dEhat_educ_p2 + dZ_Ehat_educ_p2;
        

%             % P3
%             X_educ_p3 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT*y_in.y1_educ_p1), X_educ_m3];
%             Z_educ_p3 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*Cons.P), X_educ_m3];
%             dim_Z_p3 = size(Z_educ_p3,2);
%             Ehat_educ_p3 = Trans_mat*y_in.y1_educ_p1 - X_educ_p3*beta_educ_p3;
%             MC_educ_p3 = bsxfun(@times,Ehat_educ_p3,Z_educ_p3);
% 
%             dEhat_educ_p3_beta = -X_educ_p3;
%             dEhat_educ_p3 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_X_p1+dim_X_p2), dEhat_educ_p3_beta,  sparse(obs,dim_X_p4)];
%             Z_dEhat_educ_p3 = Z_educ_p3'*dEhat_educ_p3;
% 
%             d_educ_p3 = Z_dEhat_educ_p3;
        
            % P4
            X_educ_p4 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT*y_in.y1_educ_p4), X_educ_m4];
            Z_educ_p4 = [Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*Cons.P), Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*A_starts), Trans_mat*bsxfun(@times,LMH_educ,W_ij_MAT^2*Cons.y0_educ), X_educ_m4];
            dim_Z_p4 = size(Z_educ_p4,2);

            Ehat_educ_p4 = Trans_mat*y_in.y1_educ_p4 - X_educ_p4*beta_educ_p4;
            MC_educ_p4 = bsxfun(@times,Ehat_educ_p4,Z_educ_p4);

            dEhat_educ_p4_beta = -X_educ_p4;
            dEhat_educ_p4_A = -(bsxfun(@times,Trans_mat*LMH_educ(:,1),beta_educ_p4(13)*d_A + beta_educ_p4(16)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,2),beta_educ_p4(14)*d_A + beta_educ_p4(17)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_educ(:,3),beta_educ_p4(15)*d_A + beta_educ_p4(18)*d_W_ij_A) );
            dEhat_educ_p4 = [sparse(obs,dim_BGD), dEhat_educ_p4_A, sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_X_p1+dim_X_p2+dim_X_p3), dEhat_educ_p4_beta];
            Z_dEhat_educ_p4 = Z_educ_p4'*dEhat_educ_p4;

            dZ_Ehat_educ_p4_A = [sparse(3,Cons.n); ...
                Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_W2_ij_A); Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_W2_ij_A); Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_W2_ij_A); ...
                sparse(12,Cons.n); ...
                Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_A); Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_A); Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_A); ...
                Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,1),d_W_ij_A); Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,2),d_W_ij_A); Ehat_educ_p4'*bsxfun(@times,Trans_mat*LMH_educ(:,3),d_W_ij_A); ...
                sparse(6,Cons.n)];

            dZ_Ehat_educ_p4 = [sparse(dim_Z_p4,dim_BGD), dZ_Ehat_educ_p4_A, sparse(dim_Z_p4,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];

            d_educ_p4 = Z_dEhat_educ_p4 + dZ_Ehat_educ_p4;
            
            %%% GENDER
            beta_gender = starts(dim_BGD+Cons.n+2*dim_alpha+dim_beta+1:dim_BGD+Cons.n+2*dim_alpha+dim_beta+dim_beta);
            dim_X_p1 = 0;
            dim_X_p2 = 0;
            dim_X_p3 = 0;
            %beta_gender_p3 = beta_gender(1:18);
            %dim_X_p3 = size(beta_gender_p3,1);
            beta_gender_p4 = beta_gender(1:24);
            dim_X_p4 = size(beta_gender_p4,1);              
            
            %         % P1
            %         X_gender_p1 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT*y_in.y1_gender_p1), X_gender_m1];
            %         Z_gender_p1 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*Cons.P), X_gender_m1];
            %         dim_Z_p1 = size(Z_gender_p1,2);
            %         Ehat_gender_p1= Trans_mat*y_in.y1_gender_p1 - X_gender_p1*beta_gender_p1;
            %         MC_gender_p1 = bsxfun(@times,Ehat_gender_p1,Z_gender_p1);
            %         
            %         dEhat_gender_p1_beta = -X_gender_p1;
            %         dEhat_gender_p1 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), dEhat_gender_p1_beta,  sparse(obs,dim_X_p2+dim_X_p3+dim_X_p4)];
            %         Z_dEhat_gender_p1 = Z_gender_p1'*dEhat_gender_p1;
            %         
            %         d_gender_p1 = Z_dEhat_gender_p1;


            %         % P2
            %         X_gender_p2 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT*y_in.y1_gender_p2), X_gender_m2];
            %         Z_gender_p2 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*Cons.P), Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*A_starts), X_gender_m2];
            %         dim_Z_p2 = size(Z_gender_p2,2);
            %         
            %         Ehat_gender_p2 = Trans_mat*y_in.y1_gender_p2 - X_gender_p2*beta_gender_p2;
            %         MC_gender_p2 = bsxfun(@times,Ehat_gender_p2,Z_gender_p2);
            %         
            %         dEhat_gender_p2_beta = -X_gender_p2;
            %         dEhat_gender_p2_A = -(bsxfun(@times,Trans_mat*LMH_gender(:,1),beta_gender_p2(13)*d_A + beta_gender_p2(16)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,2),beta_gender_p2(14)*d_A + beta_gender_p2(17)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,3),beta_gender_p2(15)*d_A + beta_gender_p2(18)*d_W_ij_A) );
            %         dEhat_gender_p2 = [sparse(obs,dim_BGD), dEhat_gender_p2_A, sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_X_p1), dEhat_gender_p2_beta, sparse(obs,dim_X_p3+dim_X_p4)];
            %         Z_dEhat_gender_p2 = Z_gender_p2'*dEhat_gender_p2;
            %         
            %         dZ_Ehat_gender_p2_A = [sparse(3,Cons.n); ...
            %             Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_W2_ij_A); Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_W2_ij_A); Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_W2_ij_A); ...
            %             sparse(9,Cons.n); ...
            %             Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_A); Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_A); Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_A); ...
            %             Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_W_ij_A); Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_W_ij_A); Ehat_gender_p2'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_W_ij_A)];
            %  
            %         dZ_Ehat_gender_p2 = [sparse(dim_Z_p2,dim_BGD), dZ_Ehat_gender_p2_A, sparse(dim_Z_p2,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
            %          
            %         d_gender_p2 = Z_dEhat_gender_p2 + dZ_Ehat_gender_p2;            
            
            
        
%             % P3
%             X_gender_p3 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT*y_in.y1_gender_p1), X_gender_m3];
%             Z_gender_p3 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*Cons.P), X_gender_m3];
%             dim_Z_p3 = size(Z_gender_p3,2);
%             Ehat_gender_p3 = Trans_mat*y_in.y1_gender_p1 - X_gender_p3*beta_gender_p3;
%             MC_gender_p3 = bsxfun(@times,Ehat_gender_p3,Z_gender_p3);
% 
%             dEhat_gender_p3_beta = -X_gender_p3;
%             dEhat_gender_p3 = [sparse(obs,dim_BGD), sparse(obs,Cons.n), sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_X_p1+dim_X_p2), dEhat_gender_p3_beta,  sparse(obs,dim_X_p4)];
%             Z_dEhat_gender_p3 = Z_gender_p3'*dEhat_gender_p3;
% 
%             d_gender_p3 = Z_dEhat_gender_p3;

            % P4
            X_gender_p4 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT*y_in.y1_gender_p4), X_gender_m4];
            Z_gender_p4 = [Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*Cons.P), Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*A_starts), Trans_mat*bsxfun(@times,LMH_gender,W_ij_MAT^2*Cons.y0_gender), X_gender_m4];
            dim_Z_p4 = size(Z_gender_p4,2);

            Ehat_gender_p4 = Trans_mat*y_in.y1_gender_p4 - X_gender_p4*beta_gender_p4;
            MC_gender_p4 = bsxfun(@times,Ehat_gender_p4,Z_gender_p4);

            dEhat_gender_p4_beta = -X_gender_p4;
            dEhat_gender_p4_A = -(bsxfun(@times,Trans_mat*LMH_gender(:,1),beta_gender_p4(13)*d_A + beta_gender_p4(16)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,2),beta_gender_p4(14)*d_A + beta_gender_p4(17)*d_W_ij_A) + bsxfun(@times,Trans_mat*LMH_gender(:,3),beta_gender_p4(15)*d_A + beta_gender_p4(18)*d_W_ij_A) );
            dEhat_gender_p4 = [sparse(obs,dim_BGD), dEhat_gender_p4_A, sparse(obs,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4), sparse(obs,dim_X_p1+dim_X_p2+dim_X_p3), dEhat_gender_p4_beta];
            Z_dEhat_gender_p4 = Z_gender_p4'*dEhat_gender_p4;

            dZ_Ehat_gender_p4_A = [sparse(3,Cons.n); ...
                Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_W2_ij_A); Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_W2_ij_A); Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_W2_ij_A); ...
                sparse(12,Cons.n); ...
                Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_A); Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_A); Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_A); ...
                Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,1),d_W_ij_A); Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,2),d_W_ij_A); Ehat_gender_p4'*bsxfun(@times,Trans_mat*LMH_gender(:,3),d_W_ij_A); ...
                sparse(6,Cons.n)];

            dZ_Ehat_gender_p4 = [sparse(dim_Z_p4,dim_BGD), dZ_Ehat_gender_p4_A, sparse(dim_Z_p4,dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];

            d_gender_p4 = Z_dEhat_gender_p4 + dZ_Ehat_gender_p4;
            
        else
            MC_educ_p3 = [];
            MC_educ_p4 = [];
            d_educ_p3 = [];
            d_educ_p4 = [];

            MC_gender_p3 = [];
            MC_gender_p4 = [];
            d_gender_p3 = [];
            d_gender_p4 = []; 
        end
        
        
        
        MC_educ_p1 = [];
        MC_educ_p2 = [];
        d_educ_p1 = [];
        d_educ_p2 = [];
        MC_gender_p1 = [];
        MC_gender_p2 = [];
        d_gender_p1 = [];
        d_gender_p2 = [];
        
        MC_educ_p3 = [];
        d_educ_p3 = [];
        MC_gender_p3 = [];
        d_gender_p3 = [];
        

        MC_educ = [MC_educ_m1, MC_educ_m2, MC_educ_m3, MC_educ_m4, MC_educ_p1, MC_educ_p2, MC_educ_p3, MC_educ_p4];
        Grad_educ = [d_educ_m1; d_educ_m2; d_educ_m3; d_educ_m4; d_educ_p1; d_educ_p2; d_educ_p3; d_educ_p4];
        %Grad_educ = [Grad_educ(:,1:dim_BGD+Cons.n), sparse(size(Grad_educ,1),dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta), Grad_educ(:,dim_BGD+Cons.n+1:dim_BGD+Cons.n+dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        Grad_educ = [Grad_educ, sparse(size(Grad_educ,1),dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        
        MC_gender = [MC_gender_m1, MC_gender_m2, MC_gender_m3, MC_gender_m4, MC_gender_p1, MC_gender_p2, MC_gender_p3, MC_gender_p4];
        Grad_gender = [d_gender_m1; d_gender_m2; d_gender_m3; d_gender_m4; d_gender_p1; d_gender_p2; d_gender_p3; d_gender_p4];
        %Grad_gender = [Grad_gender(:,1:dim_BGD+Cons.n+dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta), sparse(size(Grad_educ,1),dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];
        Grad_gender = [Grad_gender(:,1:dim_BGD+Cons.n), sparse(size(Grad_educ,1),dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta), Grad_gender(:,dim_BGD+Cons.n+1:dim_BGD+Cons.n+dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta)];


        %G2 = [MC_educ, MC_gender];
        G2 = [MC_educ, MC_gender];

        %Grad2 = [Grad_educ; Grad_gender];
        Grad2 = [Grad_educ; Grad_gender]/obs;
    
        W = speye(size(G1,2)+size(G2,2));
        G12 = [G1, G2];    
        G12_mean = mean(G12);
        
        %Grad1 = [Grad1, sparse(size(Grad1,1),120)];
        Grad1 = [Grad1, sparse(size(Grad1,1),2*(dim_X_m1+dim_X_m2+dim_X_m3+dim_X_m4+dim_beta))];
        Grad12 = [Grad1; Grad2];
        
        Qn = full(G12_mean*W*G12_mean');
        Grad_Qn = full(2*Grad12'*W*G12_mean');
        
        

        % Calculate Variance allowing for SEs to be clustered within
        % schools
        S = zeros(size(G12,2),size(G12,2));
        k1 = 0;
        for s=1:Cons.schools
            %s
            sch_size_s = Cons.sch_size(s);
            S = S + G12(k1+1:k1+sch_size_s*(sch_size_s-1),:)'*G12(k1+1:k1+sch_size_s*(sch_size_s-1),:)/obs^2;
            k1 = k1 + sch_size_s*(sch_size_s-1);
        end

        Z = inv(Grad12'*W*Grad12);
        %Var_GMM = inv(Grad12'*Grad12) * (Grad12'*S*Grad12) * inv(Grad12'*Grad12)';
        Var_GMM = Z * (Grad12'*W*S*W*Grad12) * Z';
        

    end


    
end

